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Abstract. The standard cold dark matter (CDM) model predicts too many and too dense 
small strnctnres. We consider an alternative model that the dark matter nndergoes two- 
body decays with cosmological lifetime r into only one type of massive danghters with non- 
relativistic recoil velocity I4. This decaying dark matter model (DDM) can snppress the 
strnctnre formation below its free-streaming scale at time scale comparable to r. Comparing 
with warm dark matter (WDM), DDM can better red nee the small strnctnres while being 
consistent with high redshfit observations. We stndy the cosmological strnctnre formation 
in DDM by performing self-consistent N-body simnlations and point ont that cosmological 
simnlations are necessary to nnderstand the DDM strnctnres especially on non-linear scales. 
We propose empirical fitting fnnetions for the DDM snppression of the mass fnnetion and the 
concentration-mass relation, which depend on the decay parameters lifetime r, recoil velocity 
14 and redshift. The fitting fnnetions lead to acenrate reconstrnotion of the the non-linear 
power transfer fnnetion of DDM to CDM in the framework of halo model. Using these resnlts, 
we set constraints on the DDM parameter space by demanding that DDM does not indnee 
larger snppression than the Lyman-a constrained WDM models. We fnrther generalize and 
constrain the DDM models to initial conditions with non-trivial mother fractions and show 
that the halo model predictions are still valid after considering a global decayed fraction. 
Finally, we point ont that the DDM is nnlikely to resolve the disagreement on clnster nnmbers 
between the Planck primary CMB prediction and the Snnyaev-Zeldovich (SZ) effect nnmber 
connt for r 
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1 Introduction 

Although there are various evidences supporting the existence of dark matter, e.g. the galaxy 
rotation curves [1-3] and bullet clusters [4], the nature of dark matter is still unknown. 
Today’s concordance cosmology model, which has shown great success in matching with the 
precision measurements of the cosmic microwave background anisotropies (CMBA) [5, 6] and 
surveys of the large scale structures (LSS) [7, 8], assumes existence of cold dark matter (CDM). 
However, observations on galactic and sub-galactic scales have shown conflicts with the CDM 
predictions. Firstly, high resolution simulations suggest far more dwarf satellites in the local 
environment than have been observed [9]. Secondly, the CDM haloes have cuspy profiles, 
which contradicts with the observed shallower profiles of dwarf galaxies [10, 11]. Thirdly, the 
largest CDM Milky Way satellites have much larger circular velocities to match with observed 
satellites [12], and fourthly the observed galaxy velocity function is significantly lower than 
the CDM prediction [13, 14]. These problems indicate that the CDM model may generate too 
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many and too dense small structures. Alternative dark matter models, with the property of 
preserving the large scale virtues of CDM while suppressing small scale structures, are worth 
considering. In this study, we assume that dark matter undergoes decays. Early discussions 
of unstable dark matter can be found in Ref. [15, 16]. The decaying dark matter (DDM) can 
suppress structure formation in general, but the detailed process still depends on the explicit 
decay models. 

One well studied scenario is the relativistic decay. This model naturally requires the 
lifetime (r) to be comparable to cosmic age (Hq^) to preserve dark matter in today’s universe. 
Consequently, haloes would experience adiabatic expansion due to the slow mass loss. With 
half of the dark matter disappeared. Ref. [17] found that the expansion can significantly lower 
the halo density concentration and prohibit star formation in dwarf satellites. The relativistic 
products might also interact with the inter-galactic medium and further alter the formation 
of first structures and the reionization history of the universe [18-20]. However, relativistic 
decays should change the expansion history of the universe and hence the growth rate of 
perturbations. The lifetime of DDM in this scenario is tightly constrained by the Integrated 
Sachs-Wolfe (ISW) effect of the CMBA, so that less than 10% of the decays could have 
occurred till now [21, 22], and the upper limit might become even restrictive from upcoming 
weak lensing surveys [23]. 

These constraints do not apply if the decay products keep most mass of their mothers 
and so remain non-relativistic. The lifetime of this kind can be very short with the decays 
completed within the radiation domination epoch. The produced daughters can act effectively 
as warm dark matter (WDM) [24], where the linear theory is adequate to understand the 
effects on perturbations. With sufficient free-streaming length and low phase-space density, 
the daughters may alleviate the CDM problems [24-26], although it is still challenging if they 
are decayed from thermally produced Weakly Interacting Mass Particles (WIMPs) [27]. Ref. 
[28] also argued that the lifetime comparable to the formation of first galaxies may offer an 
explanation of the ultra-high energy cosmic rays. 

Contrary to the short-lifetime decay, we consider non-relativistic decays with long life¬ 
time (r > Hq^). We will show that such models are completely different from the WDM 
model. To be more explicit, we consider models in which the mother particle undergoes 
two-body decay with just one type of massive daughters, where the only two possibilities are 


ddm —>■ dm + I 


( 1 . 1 ) 


and 


ddm —)■ dm + dm. 


( 1 . 2 ) 


The mother particle ddm is the decaying dark matter and the daughter particle dm is the 
stable dark matter. They are assumed to be collisionless in the structure formation. In Eq. 
(1.1) (Model A), decay also produces a relativistic daughter I, which might be a photon or 
lepton of the standard model of particle physics. We further assume I to be dark to avoid 
the constraints from the gamma ray and neutrino measurement of the galactic center and 
the diffuse background [29]. This restriction is not needed in Eq. (1.2) (Model B). In both 
situations, dm receives recoil velocity in the center-of-mass frame of ddm after production. 
With a small mass difference Am <C m between the mother and the total mass of daughters, 
the recoil velocities ( 14 ) are 

Am 
Vk = c - 


- 2 


m 


(1.3) 



and 


Vi = cd^. (1.4) 

V m 

respectively, where c is the speed of light and m is the mass of ddm. We dehne the lifetime 
and recoil velocity as the decay parameters and restrict 14 to be less than 1000 km/s for this 
study. 

The DDM effects within this parameter space automatically mix with the non-linear 
gravitational evolution making simple extension of the linear calculation from short to long 
lifetime inadequate. The induced coupled equations of structure formation have to be solved 
numerically from the hrst principle. However, previous numerical studies of this problem have 
mainly focused on the properties of isolated haloes [30, 31]. Without the cosmological envi¬ 
ronment, the growth histories of structures are oversimplihed and the statistical information 
such as the mass function of haloes and matter power spectrum cannot be determined. In 
this work, we report cosmological simulations of this problem and develop empirical relations 
to quantify the DDM effects on non-linear scales. 

We organize the paper as follows: the equations of the DDM structure formation and 
the N-body implementation are described in Section 2. In Section 3, we present the features 
of DDM on structure formation, which are originated separately from the decay parameters. 
In Section 4, we model the power spectrum suppression of DDM to CDM in the framework 
of halo model with unbounded mass. Empirical functions are proposed for the DDM mass 
function and concentration-mass (c-M) relation. In Section 5, we discuss the constraints of 
the DDM parameter space by comparing the power suppression with that of the Lyman-a 
limited WDM models. In this part, we also relax our implicit assumption that the mother 
particles completely dominate the dark matter initially. Furthermore, we argue that the 
DDM suppression is not promising as an explanation of the Planck cluster count disagreement 
[32, 33]. Finally, we summarize in Section 6. 


2 The N-body method and simulations 
2.1 Equations of the structure formation 


In the non-relativistic situation. Am/m is typically smaller than the order of 10“^. It is thus 
a good approximation to set the daughter particle’s mass to be the same as the mother’s for 
Model A and half of the mass for Model B. We use the comoving coordinate x and dehne the 
particle momentum as 

2 

p = am^ — , (2.1) 


where a is the scale factor and mx is the mass of a particle. The dynamics of the momentum 
obeys 


dp 

dt 


Vh4>, 

a 


( 2 . 2 ) 


where <5<h is the peculiar potential satisfying 


= 4ttG [/9(x, t) — p] . 


(2.3) 


A is 


The Boltzmann equation for the distribution function of the mother particle /m in model 


d/M(x,p,f) 

dt 


-A/M(x,p,t), 


(2.4) 
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See Section 3.1. 






where A = ln2/r and ^ = + + The distribution function of the daughter 

particles satishes 


d/D^(x,P,t) 

dt 


= j A/M(x,p',t)(5D(|p'- p| - aml4)^d^ 


P • 


(2.5) 


The delta function states the selection rule that the mother should have a velocity difference 
14 with its daughter. Because the decays are isotropic in the rest frame of the mother, the 
constant A is calculable from the number conservation: 


f d/M(x,p,t) 

J dt ^ 


dt 


d^p = 0, 


( 2 . 6 ) 


and we obtain 

A = d'K{amVkf‘. (2-7) 

For Model B, the Boltzmann equation of the mother particles is unchanged. But there 
are a few modihcations for that of the daughter’s. Firstly, the selection rule is altered to 

Ip' — 2p| — amVk = 0, (2.8) 


as the daughter takes only half of the mother’s mass. Besides, the number of produced massive 
daughters is also doubled in each decay. The equation is then 


dt 


= J 2A/M(x,p',t)hij(|p'- 2p| - aml4)^d3p'. (2.9) 


Similarly, the constant B is determined from a modihed number conservation relation 

p{2), 


d/M(x,p,t) l d/)j4x,p,t) 
dt 2 dt 


d'^p = 0, 


leading to 


B = 


Tr{amVk)^ 


( 2 . 10 ) 


( 2 . 11 ) 


The solutions of the daughter’s Boltzmann equations in the two models are actually 
related. It is easy to verify that the solution of Eq. (2.9) can be expressed using that of Eq. 
(2.5), if 

/B^(x,p,t) = 16/i5^^(x,2p,t). (2.12) 


The density held of Model B is then 


p(^)(x,t) = J |^m/M(x,p,t) + 

= t), 


•( 2 ), 


d^p 


(2.13) 


from which we have shown that Model A and B are actually identical in structure formation 
given the same decay parameters, although their particle physics are different. We therefore 
will not distinguish them hereafter. 
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Figure 1. An illustration of the N-body description of the DDM models. In the left, we show a 
sampling of a local density with 20 mother particles at time t. They can have very similar momenta if 
the phase-space element is defined with a narrow A^p. In the right, we show the daughter simulation 
particles after a short time of decaying. The daughters are split from the mother particles of the 
decayed mass and recoiled randomly. The relative positions of the daughters to the mothers represent 
the directions of the recoil, and we have shown only one daughter simulation particle for each mother. 
In both plots, the particles are shown in the frame of the mothers. 


2.2 N-body algorithm of DDM simulations 

N-body simulation is essentially a method to track the underlying phase-space evolution of 
the matter field. In a local region, a collection of microscopic particles with similar velocities 
defines a phase-space element which is then represented by several simulation particles. 
We illustrate this in the left panel of Fig. 1, where a phase-space element is sampled by 
20 simulation particles of similar momenta. In the CDM scenario, the evolution that let 
each simulation particle follow the dynamic equation Eq. (2.2) is equivalent to solving the 
collisionless Boltzmann equation. In the DDM scenario, microscopic daughter particles after 
decays leave the original phase-space element of the mothers with a constant recoil velocity 
to all possible directions. Therefore, sampling the new phase-space elements of the daughters 
with simulation particles is equivalent to solving the Boltzmann equations of Eq. (2.4) and 
Eq. (2.5). In the right panel of Fig. 1, we approach this by letting each simulation particle 
split the decayed mass to a new simulation particle with a randomly directed velocity. The 
splitting and new-born particles are noted as the mother and daughter simulation particles. 
The density field inferred from the daughter simulation is proportional to the density field 
inferred from the mother particles, representing that the decays of the microscopic particles are 
proportional to the local density. Also the dispersion of the microscopic daughter particles is 
represented after considering the splittings of all the mother simulation particles. Notice that 
the daughter simulation particles represent the mass emitted from the decays, they are not 
further split in the following evolution. This N-body description also explains the equivalence 

^It equals to ma;/a;(x, p, t)A®xA^p, for particles of mass rrix in a space interval A®x and a momentum 
interval A^p of the Boltzmann function A(x,p, t). 
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of Model A and B, because they generate the same isotropic daughter phase-space elements 
in the frame of the mothers. 

Previously, Peter et al. [31] proposed an N-body algorithm to study the DDM effects 
on isolated haloes, which has recently been applied to cosmological structures [34, 35]. Pe¬ 
ter’s method preserves the total number of simulation particles and samples new phase-space 
elements of daughters by randomly choosing and kicking the mother simulation particles ac¬ 
cording to the decay probability of a time interval. Their algorithm and ours shall be identical 
for the well sampled phase-space elements of mothers. However, for the badly resolved ones, 
such as the small structures or the inner regions of haloes, the variance of random picking 
from the expected value will increase quickly as the number of mother simulation particles 
drops {ojn oc l/\/]V). As a result, the lifetime is effectively not uniform in Peter’s method, 
while ours ensures that. 

To describe our algorithm, we introduce two additional parameters: the number of 
daughter simulation particles Ng produced at each split and the number of splittings fg. Ng 
means how well the velocity dispersion is sampled at a local position and fg determines the 
accuracy of updating the decayed mass. For a system evolves to time Tg, the split only 
happens when the decays accumulate to the mass fraction 

f ln2 Tg\ 

r/ = l-expf —— ■ yY (2.14) 

In actual simulations, we only use fg such that rj is just a few percent. These additional 
parameters are artificial. We check the convergence of our simulation over their choices in 
Appendix A.l and conclude that fg = 10 and A^^ = 1 are adequate for our study. 

Our N-body method can be generalized to other decay channels and their mixture with 
branching ratios. For instance, if the decay involves multi-massive daughters as 

ddm ^ dmi + dm 2 i (2.15) 

the mother simulation particle needs to split into two types of the daughter simulation par¬ 
ticles, and for three-body decays, such as 

ddm —)• dmi + dm 2 + ^ (2.16) 

simulation can be also be made after adjusting the amplitudes of the recoil velocities to follow 
certain distributions. 

We implement the method in the public N-body code Gadget2 [36]. We assign the 
daughter simulation particles a larger gravitational smoothing length to prevent the two-body 
relaxations when the lighter daughters are close to the heavy mothers. Other modifications 
are made in the tree-code to ensure the accuracy of the force with different softening lengths. 
Unique IDs are also assigned to the daughter simulation particles so that their evolution can 
be traced. 

2.3 N-body simulations 

We assume a flat universe with the cosmological parameters Vim = 0.3, Vlb = 0.049, Da = 
0.7, h = 0.7, Ug = 0.96 and cjg = 0.8, which are in between the WMAP-9yr [37] and 
the Planck 2013 results [6]. In the matter dominated epoch, we neglect the baryons and 
assume DM makes up all the mass fraction of Vim and so do the simulation particles. A 
correction of the baryon presence will be considered later. The decay parameters are set as 
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Label 

Gyr 


L 

fs 

Ng 

ri/Ng 

km/s 

Mpc//i 

DS-1 

13.48 

100 

50, 20 

10 

1 

6.70% 

DS-2 

26.20 

100 

50, 20 

10 

1 

3.50% 

DS-3 

13.48 

200 

100, 50 

10 

1 

6.70% 

DS-3a 

13.48 

200 

50 

24 

1 

2.85% 

DS-3b 

13.48 

200 

50 

10 

2 

3.35% 

DS-4 

26.20 

200 

100, 50 

10 

1 

3.50% 

DS-5 

13.48 

500 

256, 100 

10 

1 

6.70% 

DS-6 

26.20 

500 

256, 100 

10 

1 

3.50% 

DS-7 

13.48 

1000 

256 

10 

1 

6.70% 

DS-8 

26.20 

1000 

256 

10 

1 

3.50% 

WS-0.5 

- 

- 

50 

- 

- 

- 


Table 1. Columns from left to right: run label; lifetime of DDM (r); recoil velocity (14); simulation 
box size (L), split frequency (/s); number of daughters at each split (Ng) and daughter to mother 
mass ratio (rj/Ns). DS-3a and DS-3b are runs with the larger fg and Ng, respectively. Bigger boxes 
are used for higher recoil velocities in order to capture their larger suppression scales. For each box 
size, an additional CDM simulation is run without being listed here. 


r = {13.48,26.20} Gyr, which corresponds to 50% or 30% decayed fraction at z = 0, and 
14 = {100, 200, 500,1000} km/s for the DDM simulations. 

To highlight the effects of DDM, we also perform CDM and WDM simulations for 
comparisons. The initial conditions of CDM simulations are generated at Zi = 100 with 256^ 
simulation particles using the Zeldovich approximation, where we have adopted the BBKS 
formalism [38] to calculate the transfer function. We consider the simulation box sizes of 256, 
100, 50 and 20 Mpc/h with the same number of particles to explore different mass scales. 
The WDM simulation is set up from the fitting transfer function in Ref. [39] with sterile 
neutrino of 0.5 keV. Because of the long lifetime nature of the DDM models, we neglect the 
decay effects before Zi and start the DDM simulations from the same initial conditions as 
the CDM ones. The DDM simulations also occupy different box sizes according to the recoil 
velocities under the principle of using smaller box to resolve the effects of smaller 14. We will 
see in Fig. 7 and Fig. 9 that the DDM effects are consistent in the halo mass functions and 
profiles over the simulation boxes and thus the mass resolutions. Notice that with respect 
to the DDM models, there is one more hidden parameter /j that specifies the mass fraction 
of ddm in the initial dark matter component. Here, we assume that ddm contributes to all 
primordial dark matter. We leave further discussion of this parameter in Section 5.2. More 
details about the simulations are listed in Table 1. The simulations were run on the 72-core 
cluster of CUHK. With the default choice of the artificial decay parameters, we found that 
the DDM runs are on average six to eight times slower than the CDM runs, and the ratio can 
increase linearly with higher fg and Ng. 

To identify haloes, we adopt the density based halo finder AHF [40]. The halo boundary 
is defined where the enclosed mean density is 200 times larger than the background critical 
density Pcrit- We further select trusted haloes with more than 100 simulation particles in 
the CDM and WDM simulations. For the DDM simulations, this requirement is set for the 
mother simulation particles. We also calculate the power spectrum of the structures using the 

^From example, with half of the DM decayed at 2 = 0, the decayed fractions are just 8.4 x 10“"^ at Zi and 
2.4 X 10“® at the decoupling of CMB photons. 
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Figure 2. The dependence of comoving propagation distance on the start-redshift when the end- 
redshift changes. The recoil velocity here is T4 = 1000 km/s. The solid (red), dashed (green) and 
dashed-dotted (blue) lines represent the end-redshift at Ze = 0, I and 2, respectively. The recoil 
velocity here is 1000 km/s. The result can be linearly scaled to other 14- 


nearest grid assignment on 1024^ grids, where the highest k is truncated at half of the Nyquist 
frequency to avoid the aliasing effect due to the discrete Fourier transform [41]. For DDM 
models simulated with two boxes, we prefer to use the larger box simulations to calculate 
the DDM power suppression to CDM, where the truncation of long wave modes and the 
cosmic variance are less important to the power at quasi non-linear regions, which otherwise 
can cause significant underestimation of the power at /c ~ 1 /i/Mpc in simulation box of 20 
Mpc/h. In Appendix A.2, we also test the convergence of the DDM power suppression of 
DS-3 and DS-4 in box sizes of 100 and 50 Mpc/h, and find that the differences are in a few 
percent level and tend to be even smaller at high k and high redshift. 


3 Features of the DDM structure formation 


In this section, we show two unique features of the DDM models, which are from the two 
decay variables. 





0.1 1 10 

k[h/Mpc] 


Figure 3. Top: The transfer functions of DDM to CDM at z = 0. The colors (black, red, green and 
blue) denote the recoil velocities (100, 200, 500, and 1000 km/s). The solid and dashed-dotted lines 
represent different r (13.48 and 26.20 Gyr). The arrows point to the scales from Eq. (3.3) of each 
Vfe. Bottom: The power spectrum of CDM at the same redshift, where the dashed line is the linear 
power. 


3.1 Characteristic suppression scale 

Due to the production of recoiling daughters, the growth of fluctuations shall be suppressed 
under certain scale. For a daughter produced at Zg, the recoil velocity is redshifted as 

= (3,1) 

Integrating Eq. (3.1) gives us the comoving propagation distance from Zg to Zg 


Lp — 



1 1 + z ^ 

' -- dz. 

H{z)l + Zg 


(3.2) 


where H{z) is the Hubble parameter at redshift z. In Fig. 2, we plot Lp as functions of Zg 
for different Zg- With fixed Zg, Lp does not increase monotonically with larger Zg. Instead, 
competition exists between the travel time and the redshift of the peculiar velocity. The early 
produced daughters have more time to travel, but they also experience more redshift in their 
velocities, while this situation is exactly opposite for the recently produced daughters. As a 
result, a maximum value Tmax of the propagation can be reached at certain which we define 
as the free-streaming length for the observer at Zg. We can also see that the free-streaming 
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length always increases with lower becanse the danghters prodnced at all redshifts are 
allowed to travel with more time. 

Similar as Ref. [42] for WDM, the characteristic scale from which the snppression 
begins can be estimated as 

ks ^ (3.3) 

-^max 

where kg is independent of r. To quantify the relative power spectra of DDM to CDM, we 
define the transfer fnnction of DDM as 


-t-tDDM 


{k,z) 


Pppuik, z) 
PcDMik,z) 


(3.4) 


In the top panel of Fig. 3, we show that the approximation Eq. (3.3) agrees well with the 
simnlations, especially for the high recoil velocity cases. In the bottom panel of Fig. 3, we 
fnrther mark the snppression scales on the CDM power spectrnm at z = 0. The CDM power 
begins to enter the non-linear region at /cnl — 0-3 h/Mpc. Based on /cnlj we divide the DDM 
simnlations into two sets. The small snppression set (DS-1 to DS-4) only modifies the non¬ 
linear power, while the large snppression set (DS-5 to DS-8) indnces snppression already on 
linear scales. However, in both sets, the DDM snppression is more significant toward smaller 
scales, reflecting the importance of inclnding the non-linear evolntion for the DDM models. 


3.2 Time evolution 

Another featnre is that the snppression is always larger towards lower redshift as the fraction 
of decays accnmnlates. In Fig. 4, we show the tendency by comparing the redshift arranged 
snapshots of the CDM, WDM and DDM simnlations. We see that at high redshift the strnc- 
tnres are barely discernible between CDM and DDM. At lower redshfit, althongh the DDM 
simnlation still preserves the overall scheme of the filamentary strnctnres as CDM, the dense 
regions are qnite extended. In contrast, the WDM strnctnres appears to differ from the CDM’s 
mostly at high redshift. To better show the differences, we plot the DDM and WDM power 
transfer fnnctions relative to CDM in Fig. 5. As expected, we observe more snppression in the 
DDM simnlation at lower redshift, which is opposite to the regeneration of small scale powers 
in WDM [43-45]. This opposite evolntion is intrinsic to the models and shonld be nsefnl in 
distingnishing the DDM and WDM models. Besides, the high reionization redshift (z ~ 6) 
inferred from the qnasars absorption lines [46] reqnires snfhcient small scale flnctnations at 
high redshift. Being more like the CDM model in the early epoch, the long-lifetime DDM 
models conld more easily be consistent with these observations withont compromising the 
snppression at lower redshift. In contrast, reionization alone has set considerably stringent 
constraints on the mass of WDM particles [47-49]. 


4 Modelling the DDM suppression 
4.1 The halo model of DDM 

We try to reconstruct the non-linear DDM transfer fnnctions nsing halo model. The standard 
halo model assnmes all matter in the nniverse in form of haloes. To calcnlate the power 
spectrnm, the nnmber density, spatial distribntion and profiles of haloes need to be given (see 
Ref. [50] for a detailed review). However, this assnmption does not apply for models that 
can snppress the formation of small haloes, snch as WDM and DDM of this stndy, dne to the 
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logio(l + 6) 

Figure 4. Snapshots of 50 Mpc//i width and 10 Mpc/h thickness of the CDM, DDM and WDM 
simulations. Rows from top to bottom correspond to redshift at z=0, 2 and 4. The DDM simulation 
DS-3 and WDM simulation WD-0.5 are used for the plot. 

fact that there is always unbounded mass. To deal with this problem, the halo model was 
extended in Ref. [51]. We follow their method and briefly discuss the halo model of this kind. 
The idea is to separate the density field into two parts, 

Pm{x) = Ph{x) + Psix), (4.1) 

where ph and ps are the densities of halo and smooth mass. Averaging over the volume, the 
mean density is pm = Ph~\~ Ps- The halo contribution to the average density is related to its 


CDM 


DDM 


WDM 
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Figure 5. The evolution of the WDM and DDM transfer functions from simulations in Fig. 4, 
shown at ^=4, 2 and 0 with the solid (red), dashed (green) and dot-dashed (blue) lines respectively. 
The arrows represent the evolution tendencies from high to low redshift. 


mass function as 


Ph = 


roc 

J M, 


dMn{M)M, 


(4.2) 


where n{M) = dA^(> M)/dM is the halo mass function and Mcut is a cutoff mass below 
which haloes are expected not to exist. The halo mass fraction and the density contrasts of 
the two components are defined as 


and 


/ ~ Ph/Pm■, 

(4.3) 

X Px~ Px 

"X ~ - 5 

(4.4) 


Px 

where y = {h, s} stands for the halo or the smooth mass. The total density contrast is then 

5 = f6h + {l- f)Ss. (4.5) 


In the statistically homogeneous and isotropic universe, the power spectrum can be expressed 
as 

Pssik) = (1 - ffPsM + 2/(1 - f)Psh{k) + fPhhik). (4.6) 

The halo powers are decomposed into the normal one- and two-halo terms 

Phh{k) = Pih{k) + P2h{k), (4.7) 
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with the explicit form 


Pih{k) 


P2h{k) 


1 

^ / dMn{M)M‘^u^{k\M), 

P}l "/iWcut 


P^k) 


Ph 



dMMbi{M)n{M)u{k\M) 



(4.8) 


where u{k\M) is the Fourier transform of the mass normalized halo profile, bi{M) is the 
linear bias for halo of mass M, and P\in{k) is the linear power spectrum. The smooth mass 
is mapped to the underlying density field with a constant bias as <5^ ~ bgd and is assumed 
to correlate with itself and the halo density linearly. The smooth-smooth and smooth-halo 
terms are then 

Pssik) = blPunik), (4.9) 

and 

Psh{k) = r dMMbi{M)n{M)u{k\M). (4.10) 

Ph -^iVfcut 

The bias of the smooth matter is actually not a free parameter but shall be constrained by 
Eq. (4.5), since the density contrast from haloes can also be expressed as 


1 

Sh = — dMMn{M)bi{M)6, (4.11) 

Ph J M cut 

at large scales. Substituting Eq. (4.11) back in Eq. (4.5), we thus obtain 


bs = (4.12) 

where we have introduced the effective bias 

1 r°° 

bes=— dMMn{M)bi{M). (4.13) 

Ph JMcut 


To use the model, we need to understand the mass function, halo profiles and linear-bias of 
the DDM haloes at first. We examine these ingredients in the following one by one. 


4.2 The mass function 

The mass function of CDM is well studied with the excursion set theory [52], where the 
number density of haloes is related to the appearance probability of density peaks in the halo 
patch averaged over all ensembles. The usual parameterization of the mass function is 


n{M) 


IPrriff d log 0 -^ 
dlogM 


;u = 


bc{z) 

a{My 


(4.14) 


where pm is the average matter density, 5c{z) = 1.686/Il(z) is the collapse threshold and D{z) 
is the linear growth factor. The variance of the overdensity at radius R = (SM/dvrpm)”^'^^ is 
defined as 


—P^^{k)W\kR), 


(4.15) 
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Figure 6. The mass functions of the ST formalism (solid line) and measurements from CDM (round 
black points) and two DDM simulations (square and diamond points). The CDM data points are 
measured in simulation boxes of 100, 50 and 20 Mpc//i, while data points of DS-1 and DS-3 are 
obtained by combing the two simulation boxes in Table 1. The dashed lines are following the best fit 
of Eq. (4.17) with the arrows pointing to the cutoff mass indicated in Peter et al. [31] for the two 
sets of decay parameters. 


where W{y) = 3(sin?/ — ycosy)/y^ is the Fourier transform of the top-hat windows function. 
We adapt the ST formalism [53] for CDM 

/(^) = [1 + exp 

with p = 0.3, q = 0.707 and the normalization parameter A = 0.3222. 

Fig. 6 compares the ST formalism with the measured mass functions of CDM and DDM 
simulations. The ST formalism shows good agreement with the data of CDM. But for DDM, 
there is clear suppression below certain mass. DS-3 here has larger recoil velocity than DS-1, 
and so it deviates from CDM at higher mass. We can also check the predictions of previous 
isolated studies, where the haloes having the escape velocities equalling to recoil velocities 
have the mass of 8.2X10^° and 6.6x10^^ M 0 //ifor 14 = 100,200km/s, respectively. Different 
from the suggestion of Peter et al. [31] that haloes with escape velocity smaller than 14 should 
be destroyed, we observe no truncation of the DDM halo mass function below these scales 
in cosmological simulations. Oversimplihcation of the formation history in isolated studies 
shall be the reason of the difference. Since small haloes in cosmological simulation are formed 


14 - 







- 2-1012 -10123 


logio[M/Mf(Vk,z)] logio[M/Mf(Vk,z)] 

Figure 7. The mass function ratio of DDM to CDM as a function of the Mf normalized halo mass 
for different DDM parameters and redshifts. The dash lines represent the best fit of Eq. (4.17). The 
filled points refer to the measurements in larger boxes, and the not-filled points are from the smaller 
boxes as in Table 1 for DS-1 to DS-4, where consistency of simulation boxes and resolutions can be 
observed. 


earlier, they can survive the decays through adiabatic mass loss but without being completely 
destructed. Their mergers would still hierarchically form bigger haloes. 

To describe the mass function of our simulations, we develop a htting function in the 


form 


nDDM(M,z) _ ^ 

ncDM(M,z) V 


(4.17) 


in which am and Pm are htting parameters. Here, we have introduced two effective variables 


dvr 


~ 2 PmLmuxi^k, z) 


(4.18) 


and 


/d = 1 - exp 


in 2 
r 


T{z) , 


(4.19) 


where Mf is the characteristic mass of the free-streaming length Lmax and fd is the decayed 
fraction at time T of redshift z. They are designed to separate the dependence of the sup¬ 
pression on the decay parameters. Our htting function does not have explicit dependence on 
the redshift, whose inhuence is already embedded in Mj and fd- 

Fig. 7 shows the mass function ratio of DDM to CDM versus the Mf normalized halo 
mass. The compared DDM simulations in each panel have the same fd but different Mf. The 
dashed lines are the best ht of Eq. (4.17) with the parameters: am = 0.526 and Pm = 7.61. 
The htting function shows good agreement with simulations of different combinations of decay 
parameters as well as redshift. Here, we have only considered the small suppression set. The 
large suppression set is found to deviate signihcantly from this result by affecting already the 
linear power. Currently, there is no theoretical work on the mass function of DDM models. 
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Figure 8. The dwarf halo profiles randomly selected from simulations DS-1 to DS-4 at 2 ; = 0. The 
solid (red) lines are the best-fit of NFW and the vertical dash line marks the spatial resolution of each 
simulation. 


Given the simplicity of this empirical approach, our result shall motivate further theoretical 
considerations on the mass function of DDM models, for example a both scale- and time- 
dependent barrier in the excursive set theory. 


4.3 The halo profile 

The density distribution of a CDM halo is usually modelled by the Navarro-Frenk-White 
(NFW) profile [54, 55] 

12 > (4-20) 

r(l -h r/rsY 

where ps and are the normalization density and characteristic length scale. It is usually 
convenient to reparameterize the profile by two other parameters, the virial mass M and 
the halo concentration c = rvirjrs, where the virial radius r^r is related to the virial mass 
through M = 47r/3r[]j^A/9crit with A = 200 for our halo definition. These two parameters are 
related after integrating the profile to virial radius, which gives the relation 


M = A-KPsV^s 


ln(l -k c) 


c 

1 -k C 


(4.21) 


In DDM models, haloes are kinetically heated from the recoiling daughters. Previous 
semi-analytic study expected the formation of cores in dwarf haloes [56]. In Fig. 8, we ran¬ 
domly select the dwarf mass range haloes in simulations DS-1 to DS-4 to test this prediction. 


- 16 - 
















- 2-1012 -10123 


logio[M/Mf(Vk,z)] 


logio[M/Mf(Vk,z)] 


Figure 9. Ratio of the NFW concentration between DDM and CDM as a function of normalized 
mass. The dash lines are the best-fit of Eq. (4.22). The point styles are the same as in Fig. 7. 
Consistency of the DDM halo concentration can also be observed over the simulation boxes. 


We show the best-fit NFW profiles with the red lines, and nsed the vertical dashed lines to 
present the spatial resolntions of the simnlations. We observe no evidence of core developing 
within onr resolntion limits. Obvionsly, to determine whether there are cores or how they 
depend on the decay parameters is beyond the scope of these simnlations. We leave these 
qnestions to fntnre high resolntion stndies. In this paper, we still stick to the NFW profile 
for DDM haloes; however, we examine how the concentration changes with the DDM param¬ 
eters as a qnantification of the halo inner density decrease. Similar to the mass fnnction, we 
parameterize the DDM concentration as 


cddm(M,z) 

ccdm(M,z) V ) 


(4.22) 


In Fig. 9, we show the best-fit fnnction with Uc = 0.271 and 13c = 20.4. The decays make 
little change to the massive haloes (M ^ Mf) becanse of their deep gravitational potentials. 
While in the low-mass limit, the DDM concentration-mass relation (c-M relation) scales as 
c oc , indicating a tnrn-over of the concentration with significant decays (i.e. 

fd > 0.22). We also investigated the variance of the logarithmic concentration of DDM 
haloes and fonnd no obvions difference from the CDM case. 


4.4 The linear halo bias 


The linear bias of CDM haloes is nsnally understood from the peak-backgronnd split [57]. 
For the ST formalism, the bias is [53] 


ftcDM — 1 + 


- 1 


-b 


2p 


Sciz) dc{z)[l + {qiz'^y]' 


(4.23) 


where zz depends on M and the parameters p and q are the same as in Eq. (4.16). Becanse 
the DDM snppression is only important in late times, the statistics of the early density field 
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Figure 10. Left: The halo mass mapping for decay parameters r = 13.48 Gyr and 14 = 100, 200 
km/s at z = 0. Right: The ratio of the halo bias between DDM and CDM for the same decay 
parameters. The dotted lines in both panels represent the limit of infinitely small DDM suppression, 
where Mddm = -Mcdm- 


is the same as CDM. However, due to the decays, a proto-region that ought to collapse to a 
CDM halo of mass Mcdm may now form a DDM halo of mass Mddm- The two haloes could 
occupy similar positions, which can also be seen in Fig. 4. The DDM halo bias is then known 
after adapting a mass mapping between the CDM and DDM haloes. 

Physically, Mddm is expected to be smaller than Mcdm through several ways of mass 
loss: (1) recoiling daughters can directly escape from shallow potentials; (2) the decays can 
reduce halo concentrations, which makes them more vulnerable to tidal stripping; (3) with 
mean halo density reduced, the virial radius would also shrink to enclose less mass. Making an 
optimistic assumption that the decay effects still preserve the number of haloes, we construct 
the mass mapping as 

dMpDM _ racDM(McDM, z) 
dMcDM ’^DDmCMddM, 2 ) ’ 

where Mddm is a function of Mcdm and on the right hand side are the CDM and DDM mass 
functions. Because decays only modify the structures within non-linear scales for small 14, 
the linear bias of DDM halo is therefore 

6ddm(Mddm) = ^cdm(Mcdm)- (4.25) 

We solve the equations numerically by involving Eq. (4.17). In Fig. 10, we present the 
mass mapping and the bias for two sets of decay parameters. In the mass mapping, the DDM 
suppression is always stronger towards lower mass haloes. But this is not true for the bias as 
shown in the right panel, since the CDM bias is no longer sensitive to the mass difference in 
the low mass end. Instead, the largest difference appears in the medium mass range, where 
the DDM suppression is important and the CDM bias increases quickly with the halo mass. 
However, the overall difference is still small. As we will see in Fig. 11, it hardly contributes 
to the power suppression. 
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Figure 11. Ratios of the DDM power without the changes in the mass function (red solid), c- 
M relation (green dashed) and halo bias (blue dotted) to that including all the effects at z = 0. 
The decay parameters here are r = 13.48 Gyr and 14 = 200 km/s. The arrow here points to the 
characteristic suppression scale kg as defined in Eq. (3.3). 

4.5 Reconstruction of the DDM power suppression 

Before applying the halo model we described, we also need to know the mass fraction of all 
haloes. For WDM, Ref. [51] suggested a physical cutoff mass M^ut according to the WDM 
free-streaming scale. This assumption can avoid the problems of ambiguous extrapolation of 
the mass function to the low mass end and also the numerical instabilities of Eq. (4.15) when 
-^cut 0. However, we have found no sign of physical cutoff in the DDM mass function. 
However, as we show in Appendix B, the small enough haloes could act exactly like the 
smooth component. We therefore argue that the cutoff mass is still appropriate for the DDM 
halo model calculation, but now it should refer to the mass scale below which the lower mass 
haloes and the real smooth component are indistinguishable. In this sense, the standard halo 
model and the halo model with smooth component are unified. 

In Fig. 11, we firstly examine the roles of the DDM effects on the mass function, halo 
profile and halo bias on the final halo model power. Each line represents a result without 
certain modification of the halo model ingredient. We can see that the DDM reduction on the 
mass function is mostly important for the power suppression and also marks the characteristic 
scale. The change in the c-M relation starts to be important in smaller scales, while the effect 
of linear bias difference is negligible. 

In Fig. 12, we finally compare the calculated DDM transfer functions with the data from 
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Figure 12. Comparison of the modelled DDM transfer functions (lines) and the simulation data 
(points). The top and bottom panels refer to the results at z = 1 and 0. The DS-3 and DS-4 
simulations of 100 Mpc//i together with the DS-1 and DS-2 simulations of 50 Mpc/h are used for the 
plot. 


simulations. The halo model predictions describe well the data to the non-linear scales. The 
modelling seems even better if the DDM suppression is smaller. The largest mismatch is in 
simulation DS-3 at z = 0 and the range /c ~ (2 — 10) /i/Mpc with the relative error about 2%. 
We also found that the suppression tails of different lifetimes do not overlap when the scale 
is normalized by kg, meaning that the transfer function cannot be described by any function 
in the form f{k/ks,fd)- The reason is that the contributions from 14 and r are no longer 
separable for the suppression tails. We also want to draw the attention that the ingredients 
of the DDM halo model are summarized from the small suppression set simulations. A safe 
parameter range to utilize these results is 14 smaller than 200 km/s and r lager than or at 
least the same order as the cosmic age. The large suppression set simulations are found to 
have big deviations from the semi-analytical predictions. But fortunately, they are also not 
consistent with observations. 

5 Discussions 

Based on the previous results, we discuss (1) the constraints of the DDM parameter space, 

(2) the more generalized DDM models with non-trivial mother particle initial fractions and 

(3) whether the DDM suppression can be a solution to the Planck Sunyaev-Zeldovich and 
primary CMB disagreement. 
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Figure 13. The transfer functions of three DDM models compared with a 2 keV WDM at z = 3. 
The lines correspond to the DDM models of t = 10 Gyr and 14 = (50,100,200) km/s, which will 
cause less flux power of Lyman-a if the suppression is within the colored region corresponding to that 
of WDM. 


5.1 The Lyman-a constraints 

We consider the Lyman-a constraints on the DDM parameters. The Lyman-a forest is 
the neutral hydrogen absorption lines in the spectra of distant quasars (QSOs). As the 
hydrogen clouds are tracing the density perturbations, the flux spectrum Ppik, z) is imprinted 
with the information of the underlying density field at medium redshift z ~ (2 — 4) and on 
scales k ~ (0.1 — 10)/i/Mpc that have not been fully contaminated by non-linear evolution. 
Usually, the flux bias function hp{k,z) = Ppik, z)/P{k, z) that relates the flux power to 
the real density power has to be understood before interpreting the data. However, it has 
a complicated dependence on the cosmological parameters, the initial power spectrum as 
well as the parameters of the baryon physics [58-60]. A large number of hydrodynamical 
simulations are in principle needed to do so. Here we simply assume that the flux bias 
function is unchanged for the WDM and DDM models at the redshift and scales relevant to 
the Lyman-a observations. We can then translate the Lyman-a limits on WDM to DDM. 
Considering the their opposite evolution tendencies, this assumption might underestimate the 
true flux power in the DDM models and make the DDM constraints conservative. 

The preferred DDM parameters are those that cause less suppression on the CDM power 
than that from the lower mass limit of WDM particles. The Lyman-a forest has been shown 
to be sensitive to the WDM mass with the lower mass limit between (2 — 4) keV for thermal 
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Figure 14. The constraints on the DDM parameter space by translating the Lyman-a constraints 
on thermal WDM. The DDM halo model is applied in the range r > 5 Gyr and T4 < 200 km/s. The 
allowed parameter space include region A for the 4 keV WDM case, and will extend to include region 
B for the 2 keV case. Region C is ruled out for our consideration and region D is not explored as kg 
there may touch linear scales. 


WDM particles [61-64], where the data set probing higher redshift usually concludes in higher 
mass. In Fig. 13, we give an example of the translation by showing the power suppression for 
a 2 keV WDM and several DDM models at z = 3. The transfer functions of the DDM models 
are calculated from the DDM halo model. The suppression of the WDM is denoted as the 
color region, where we have adopted a well calibrated fitting formula for the WDM transfer 
function from Ref. [65]. The DDM models within the color region have larger suppression 
than the WDM model and are thus not favoured. We then survey the parameter space with 
r > 5 Gyr and I 4 < 1000 km/s. The results are shown in Fig. 14, where the parameter space 
is divided into four regions (A, B, C, and D). Region A is conservative and allowed by the 4 
keV WDM limit. Region B is in between the 4 keV and 2 keV limits. Region C is ruled out, 
but region D is still uncertain because the accuracy of the DDM halo model is inadequate. 
At r = 13.48 Gyr, the constraints of the recoil velocity are that I 4 should be smaller than 
(35 — 105) km/s. While at r = 10 Gyr, the constraints are I 4 less than (31 — 88) km/s. This 
result is consistent with a recent direct fitting of DDM models to the Lyman-a data in Ref. 
[34], where they have concluded that for r < 10 Gyr, I 4 is smaller than (30 — 70) km/s. 
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Figure 15. Comparison of the halo model predictions by replacing with fg and the transfer 
functions measured from the DDM simulations with fi =0, 0.4 and 0.8. The decay parameters here 
are t = 13.48 Gyr and 14 = 100 km/s. 


5.2 More generalized initial condition 

In all previous studies of DDM, the mother particles ddm are assumed to make up all the 
matter component initially. We relax this limit by considering that the ddm only contributes 
to parts of the initial mass, leaving the other mass as stable matter with the fraction 


f^ = l- 




(5.1) 


at Zi- The stable matter might be the daughter particles dm or baryons or even other type 
of CDM. The first possibility may be realized if ddm and dm are both WIMPs and are 
thermally produced with comparable annihilation cross sections. We expect the stable and 
unstable components are uniformly mixed initially. Effectively, if a global decayed fraction 
fg = {i — fi)fd can be mimicked by an other lifetime Tefr with /* = 0 such that 

fd{reS,z) = fg{T,z), (5.2) 

the structure formation of the two cases should be exactly the same. However, the solution 
of Eq. (5.2) is time dependent in general, unless in the very long lifetime limit (r ^ 
the effective lifetime approaches a constant 




1 -/*’ 


(5.3) 
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Figure 16. The generalized constraints of DDM models on global decayed fraction and recoil velocity. 
The regions A to D represent the same regions as those in Fig. 14. The region E bounded with the 
magenta dashed lines corresponds to the parameter space of resolving the Planck disagreement on the 
cluster number with fi = 0. 


implying that fi and r can be degenerate. 

For the normal sitnation with r ~ non-zero fi will canse a different accnmnlation 

of the decay prodnced danghters, which can be stndied by changing the r] of the N-body 
algorithm to 


^(Ti) 


f ln2 TAl (1 - fi) 

V fs) _ (1 -/i)-h ^exp (^ • Ti) 


(5.4) 


We redo the simnlation DS-1 of 50 h/Mpc with the new ratio and consider fi = 40% and 
80% that correspond to fg = 30% and 10% at ^ = 0. Since it is the new born danghters 
that canse the snppression, we also generalize the DDM halo model by replacing fd with fg. 
The transfer fnnctions from the new simulations and new halo model are compared in Fig. 
15. The good agreement shows that the global decayed fraction fg is the right parameter to 
describe the generalized DDM models. Also notice that we previonsly only snmmarized the 
DDM halo model ingredients from the simnlations with the decayed fraction of 30% and 50% 
at z = 0] the match to the new decayed fraction of 10% also demonstrates the accnracy of 
halo model to the small snppression end. Another point to notice is that to have the same 
fg at a fixed redshift, there is a degeneracy in fi and r. However, the degeneracy shonld be 
broken if the snppression is examined at other redshifts. 
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Figure 17. Halo number density differences of DDM normalized by the number density of CDM 
measured in simulations with A = 500 at z = 0. The grey region shows typical cluster mass range. 
The yellow region represents the suppression of DDM needed to resolve the Planck disagreement. 
DDM simulations that have different box sizes in Table 1 are combined to make the plot. 


Based on Fig. 14, we constrain the more generalized DDM models in Fig. 16. Fig. 16 
also allows us to make a first order correction to the DDM constraints with baryonic matter, 
which has been neglected previously. By assuming that baryons honestly trace the undecayed 
mass, their presence is equivalent to an effective fi = = 16.3%. The DDM constraints 

with all initial dark matter as mother particles at r = 13.48 Gyr are broadened to 14 less 
than (37— 120) km/s. 

5.3 Resolving the Planck disagreement? 

The Planck 2013 results have reported the constraints of the cosmological parameters using 
the number counts of the Sunyaev-Zeldovich effect selected clusters [32]. The interpolated 
parameters ag and are found to differ from those derived from the primary CMB temper¬ 
ature anisotropies, which is also indicated by the latest data release [33]. With a reasonable 
bias of the measured cluster mass to the real mass, it could lead to the number of predicted 
clusters from CMB two times larger than the observed. 

DDM may provide a natural explanation for this disagreement, as it can reduce the 
cluster number in the late evolution without interfering the CMB. Such possibility was firstly 
proposed in Ref. [66], but only linear perturbations have been done. With our N-body 
simulations, we revisit the DDM suppression on the mass function with a higher A = 500 
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for the cluster convention [32], Focusing on the cluster mass range, we see in Fig. 17 that 
the small suppression set (DS-1 to DS-4) makes no difference to CDM. The preferred DDM 
suppression then should reduce the cluster number by more than those of DS-5 and DS-6 but 
less than those of DS-7 and DS-8. For r ~ this parameter space is shown as region E 

in Fig. 16, which however is deeply inside the Lyman-a ruled-out region C. The reason is 
that to have such a large decrease of high mass haloes, smaller haloes are already suppressed 
much more. The fi here is zero. Higher /j would not change the conclusion, since that will 
require even faster decays and move the region E upwards. Also notice that our simulations 
do not have the same cosmology as Planck measured. However, the Planck cosmology will 
produce even more high mass clusters. We expect that the region E is also conservative for 
the exact Planck cosmological parameters. 

The failure of the DDM models suggests that the late-time suppression shall only occur 
in high mass objects. In this sense, mechanisms like the AGN feedback [67, 68] might be 
more plausible to resolve the disagreement. 

6 Summary and conclusion 

In this paper, we studied the cosmological structure evolution in the non-relativistic and 
long-lifetime DDM models. The decay mechanism brings unique features to the structure 
evolution, with the recoil velocity 14 determining a characteristic suppression scale and the 
lifetime r regulating the time when the suppression is important. Intrinsically different from 
the WDM, the structures in DDM models are more CDM like in early times. We argue that 
DDM models could more easily cause suppression in the late universe while being consistent 
with the high redshift observations, such as the reionization and Lyman-a forest. 

In particular, we considered the DDM models with two-body decay and only one type 
of massive daughter. The two possible cases (Model A and B) were shown to be identical 
in the structure formation. Using N-body simulations, we solved the coupled equations that 
govern the DDM structure evolution from the first principle. These cosmological simulations 
are needed to understand the effects of DDM, especially on non-linear scales. In the ana¬ 
lytical aspect, we proposed empirical functions parameterized by the characteristic mass Mf 
and decayed fraction /^, which are again functions of the decay parameters and redshift, to 
describe the DDM suppression on the mass function and the halo concentration. This also 
leads to accurate reconstruction of the non-linear power transfer function of DDM to CDM 
in the framework of halo model. The consequence of these efforts is that using the analyti¬ 
cal predictions the decay parameter space can be explored far more than the few simulated 
points. 

Additionally, we translated the constraints of WDM mass from Lyman-a forest to those 
of the DDM parameters. The results were shown in Fig. 14. For decay models with r ~ 
we found the 14 should be smaller than 105 km/s or 35 km/s for more conservative 
consideration. The DDM models are also generalized to arbitrary initial fractions of the 
mother particles. We found that the halo model is still valid after replacing with the 
global decayed fraction fg = (1 — fi)fd- Constraints were also made for the generalized DDM 
models as shown in Fig. 16. Using the constraints, we also demonstrated that the DDM 
models are unlikely to resolve the disagreement on the cluster numbers from the Planck SZ 
survey and primary CMB prediction without violating the Lyman-a limits. 

Through the study, we have shown that DDM is rich in phenomenology and its ability 
of reducing high density and suppressing the formation of small structures has been revealed. 
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Figure 18. Consistency tests of the artificial decay parameters, where DS-3a and DS-3b have either 
high value in fg and Ng than DS-3 of the same simulation box. Left side: Their mass functions at 
z = 0. The colored region represents the errors of DS-3. Right side: Their power spectra at z = 0 
with the dashed line showing the linear power. In both sides, the green and blue lines in the bottom 
panels represent the differences of DS-3a and DS-3b with respect to DS-3. 


In future work, we will substitute DDM for CDM to other dark matter related observations, 
such as weak lensing, HI surveys of the galaxy velocity function and dark matter detections. 
We also expect higher resolution simulations of our algorithm to explore more details of DDM 
structures on smaller scales Together, we may have an answer of whether DDM can be a 
better alternative of CDM. 

A Convergence tests of the simulations 
A.l Test the artificial decay parameters 

Physical properties extracted from the DDM simulation should depend on the decay param¬ 
eters (r and 14) rather than the artificial simulation parameters (/^ and Ng). As listed in 
Table 1, two test simulations DS-3a and DS-3b are performed, which have either higher fg 
or Ng than DS-3. We examine the convergence of the simulations by comparing their mass 
functions and power spectra at z = 0. As shown in the left panel of Fig. 18, DS-3a and DS-3b 
have the same number density of haloes as DS-3 in the high mass end. Meanwhile, DS-3b 
is closer to DS-3 than DS-3a in the low mass end, indicating that the parameter Ng is more 
easily converged than fg. However, the overall differences are still less than 15%. The right 
panel shows an better consistency in the power spectra, where the largest difference is less 
than 2% even on highly non-linear scales. Since changing Ng and fg will let the simulations 
run with completely different particles, we believe that the choice of fg = 10 and Ng = 1 al¬ 
ready makes the simulations converged, and we use them as the default parameters for other 
DDM simulations. 

^It would be interesting to compare the profiles and mass function of subhaloes with the zoom-in simulations 
of Wang et al. [35], where they have applied the Peter’s algorithm. 
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Figure 19. Consistency tests of the DDM transfer functions with the simulation boxes. The relative 
differences of the transfer functions from boxes of 50 Mpc/h are compared with the results from boxes 
of 100 Mpc/h. DDM parameters of DS-3 and DS-4 are shown with the red and green lines, and 
redshifts at 2 = 0 and 1 are further represented with the solid and dashed styles. 


A.2 Test the DDM suppression with box sizes 

In Fig. 19, we test the convergence of the DDM power transfer fnnctions defined in Eq. (3.4) 
of DS-3 and DS-4 with the simnlation box sizes. The resnlts of 50 Mpc//i are compared to 
the resnlts of 100 Mpc//i. The largest differences occnr in DS-3 at z = 0, which is abont 
3.5% at fc ~ 17 h/Mpc before getting smaller at higher k. The differences are also redshift 
and DDM snppression dependent, with better consistency at high redshift and with smaller 
DDM snppression. 


B Equivalence of small haloes and smooth component in the halo model 

The standard halo model can be expressed as 


P{k) = P2}iik) + Pm{k), 
where the two- and one-halo terms are 


P2Rik) = -zjrPlinik) 
Pm 


dMMbi{M)n{M)u{k\M) 


and 


0 

/*oo 


Pm{k) = ^ / dMn{M)M^u\k\M). 
Pm J 0 


The explicit form of the mass filter of a NFW halo is 

u{k\M) {sin(fcr^) [Si[(l -h c)krs] - Si(fer^)] - 


L —V--— v-'^yj (^l c)krs 

+ cos{krs) [Ci[(l -h c)krs] - Ci(A:rs)]}, 


(B.l) 

(B.2) 

(B.3) 


(B.4) 
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Figure 20. The standard halo model power calculated using different cutoff mass. In the top panel, 
the dashed, dotted, dashed-dotted and the dashed-double-dotted lines correspond to the terms in Eq. 
(B.12) and (B.13). The solid lines are the sum of the terms. The results from Mcut = 10®Mq//i 
and 1Mq//i are represented with the red and black colors, respectively. The green dotted line shows 
the upper limit of the white noise of Mcut = lO®M 0 /h. In the bottom panel, we plot the fractional 
difference of the total power of Mcut = 10®Mq//i with respect to that of Mcut = 1Mq//i. 


where Si(x) = sin(t)/tdt and Ci(x) = — cos{t)/tdt. Eq. (B.4) can reach the limit 


lim u{k\M) 

fcrs—tO 


d-irpsT^ 

M 


ln(l + c) 


c 

1 + c 


= 1 , 


(B.5) 


when krg 1, where the last equality is the mass conservation equation of the NEW prohle. 
The limit can be approached on large enough scales or in haloes of small enough r^. Demand¬ 
ing that the two-halo term equals the linear power on large scales, Eq. (B.2) and Eq. (B.5) 
lead to a non-trivial bias mass relation [69] 



dM Mhi{M)n{M) = pm- 


(B.6) 


Assuming a cutoff mass Mcut 


below which the largest wave number of interest satishes 
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^max^s (Mcut) <C 1 we have 


(•Mcut 


poo 

/ dMMbi{M)n{M)u{k\M) = I dMMn{M)bi{M) 

Jo Jo 

poo 

+ / dMMn{M)bi{M)u{k\M). 

J ^cut 

After defining the mass fraction for haloes larger than -Mcut as 

1 r°° 

f = — dMMn{M) 

Pm J Mcut 

= EIl 

— ? 

Pm 

the first term on the right hand side of Eq. (B.7) can be rewritten as 

rMcut 


>0 


dMMn{M)bi{M) = 6,(1 - /)pu 


Combining Eq. (B.9) and Eq. (B.6), we have 

1 


6 ., = 


(1 - f)Pm 

1 - 


Pm- dMMbi{M)n{M) 

Mcut 


1 -/ ’ 

with the form of the effective bias as 


b^^ = — 

Ph J Mcut 

Inserting Eq. (B.7) and Eq. (B.9) back to Eq. (B.2), we have the two-halo term 

L p poo 

P 2 R =(1 - + 2(1 - f)f^^ / dMMniM)biiM)uik\M) 

Ph JMcut 


poo 

/ dMMn{M)bi{M). 

J Mcut 


1 2 


dMMn{M)bi{M)u{k\M) 

Ph Mcut 

= (1 - /)2p„ + 2(1 - f)fPsh + fP2h. 


(B.7) 


(B.8) 


(B.9) 


(B.IO) 


(B.ll) 


(B.12) 


Using the same strategy, the one-halo term is separable as 

1 coo 1 rMcut 

Pm = f^ dMn(M)M‘^u^(k\M) + ^ dMn{M)M^ 

Ph J Mcut Pm Jo (B.13) 

= fPlh + Pn,. 


®The DDM models can also satisfy the inequality, because oc ruirjc oc 0 . 271 /d extending Eq. 

(4.22) to the low mass end. 
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Therefore, we have shown that the standard halo model and the halo model with unbounded 
mass are equivalent in calculation, except for a negligible white-noise term 


1 

7r2 

rMcut 

/ dMn{M)M‘^ 

Pm 

'o 

■ 1 

fMcut 

/ dMn(M)M 

_Pm 

Jo 

(1- 

f)^ 

pm 

A/cut 

Pm 



Mcnt 

Pm 


(B.14) 


The white noise can be smaller than 10“® (Mpc//i)^ if Mcut < lO^MQ/h. We present a 
test in Fig. 20 by calculating the standard halo model with different cutoff mass. The good 
consistency shows that haloes smaller 10^Mq//i are already indistinguishable from the smooth 
component for femax <100 h/Mpc. 
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